cx = [-1, -0.5, -0.4, -0.3, 1, 1, -1];
cy = [1, 1, 0.8, 1, 1, -1, -1];

% cx = cx + 1;
% cy = cy + 1;

[x, y] = generateBoundaryNode(cx, cy, 30);

domain2D = createDomain(x, y);

dim = size(domain2D.x, 2);

b = assemblyBoundaryVector(domain2D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

res = 30;
scale = 0.9;

[px1, py1] = meshgrid(0:res/3, 0:res/3);
[px2, py2] = meshgrid(0:res, 0:res);
[px3, py3] = meshgrid(0:res/3, 0:res/3);

px1 = reshape(px1, (res/3+1)*(res/3+1), 1) / res * 3 * 0.4 - 0.4 / 2;
px2 = reshape(px2, (res+1)*(res+1), 1) / res * 1.2 - 1.2 / 2;
px3 = reshape(px3, (res/3+1)*(res/3+1), 1) / res * 3 * 0.3 - 0.3 / 2;
py1 = reshape(py1, (res/3+1)*(res/3+1), 1) / res * 3 * 1.8 - 1.8 / 2;
py2 = reshape(py2, (res+1)*(res+1), 1) / res * 1.8 - 1.8 / 2;
py3 = reshape(py3, (res/3+1)*(res/3+1), 1) / res * 3 * 1.7 - 1.7 / 2;

shift = 1;

px1 = px1 - 0.75;
py1 = py1;
px2 = px2 + 0.35;
py2 = py2;
px3 = px3 - 0.4;
py3 = py3 - 0.1;

rx = [px1; px2; px3];
ry = [py1; py2; py3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[A, Z] = assemblyMatrix(domain2D);

theta = A \ b;

u = evaluateInteriorPoint(rx, ry, domain2D, theta);

F = evaluateInteriorDeformationGradient(rx, ry, domain2D, theta);

plotStrain( domain2D, rx, ry, F )

% plotDeformation(domain2D, rx, ry, b(1:2:2*dim), b(2:2:2*dim), u(:,1), u(:,2));